Provide in one place some code to simulate data and fit spatial capture-recapture models, for both closed and open populations.
The reference book for spatial capture-capture (SCR) models is .
Slides of an introduction to SCR models by Andy Royle,
We will go the Bayesian way. Richard Chandler provides nice slides for Bayesian inference in closed population SCR models, and open population SCR models.
We will use Nimble to fit SCR models. We might implement a few tricks to improve and reach convergence faster.
Load Nimble.
library(nimble)
We use some code shared by Jose Jimenez on the SCR Google group and adapted by Perry de Valpine.
The traps.
tr <- seq(15, 85, length = 10)
traps <- data.frame(X = rep(tr, each = length(tr)),
Y = rep(tr, times = length(tr))) # 100 coord. traps
Visualize.
viz_traps <- traps %>%
ggplot(aes(x = X, y = Y)) +
geom_point(pch = 3) +
xlim(0, 100) +
ylim(0, 100)
viz_traps
Generate population.
set.seed(10)
xlim <- c(0, 100)
ylim <- c(0, 100) # area 100 * 100 = 1e4
A <- (xlim[2] - xlim[1]) * (ylim[2] - ylim[1])/10000
A
## [1] 1
mu <- 50 # density
N <- rpois(1, mu*A) # generate population
N
## [1] 50
Generate activity centers.
s <- data.frame(s.x = runif(N, xlim[1], xlim[2]),
s.y = runif(N, ylim[1], ylim[2]))
Visualize.
viz_traps_ac <- viz_traps +
geom_point(data = s, aes(x = s.x, y = s.y), pch = 16, color = "red")
viz_traps_ac
Generate detections.
sigma <- 5
lambda0 <- 0.4
J <- nrow(traps) # nb of traps
K <- 5 # nb capture occasions
yy <- array(NA, c(N, J, K))
for(j in 1:J) {
dist <- sqrt((traps$X[j] - s$s.x)^2 + (traps$Y[j] - s$s.y)^2)
lambda <- lambda0 * exp(-dist^2 / (2 * sigma^2))
for(k in 1:K) {
yy[,j,k] <- rpois(N, lambda)
}
}
n <- apply(yy, c(2,3), sum)
Plot detections.
tot <- apply(n, 1, sum)
dat <- data.frame(traps, tot = tot)
viz_traps_ac +
geom_point(data = dat, aes(x = X, y = Y, size = tot), alpha = 0.3) +
scale_size(range = c(0, 20)) +
labs(x = "",
y = "",
size = "# detections")
Define the model.
code <- nimbleCode({
sigma ~ dunif(0, 10)
lam0 ~ dunif(0, 5)
psi ~ dbeta(1, 1)
for(i in 1:M) {
z[i] ~ dbern(psi)
s[i,1] ~ dunif(xlim[1], xlim[2])
s[i,2] ~ dunif(ylim[1], ylim[2])
dist[i,1:J] <- (s[i,1] - X[1:J,1])^2 + (s[i,2] - X[1:J,2])^2
lam[i,1:J] <- exp(-dist[i,1:J] / (2 * sigma^2)) * z[i]
}
for(j in 1:J){
bigLambda[j] <- lam0 * sum(lam[1:M,j])
for(k in 1:K) {
n[j,k] ~ dpois(bigLambda[j])
}
}
N <- sum(z[1:M])
})
Define constants, data and inits.
M <- 200
constants <- list(M = M,
K = K,
J = J)
n1 <- apply(n, 1, sum)
data <- list(n = n,
X = traps,
xlim = xlim,
ylim = ylim)
s <- cbind(runif(M, xlim[1], xlim[2]),
runif(M, ylim[1], ylim[2]))
z <- rep(1, M)
inits <- list(sigma = 0.5,
lam0 = 0.1,
s = s,
z = z,
psi = 0.5)
Build R model (not compiled yet).
Rmodel <- nimbleModel(code = code,
constants = constants,
data = data,
inits = inits)
Check whether the model is fully initialized. If you failed at providing initial values for some parameters (e.g. \(\psi\)), you’ll get NAs.
Rmodel$calculate()
## [1] -8089.009
Now compile the model in C++.
Cmodel <- compileNimble(Rmodel)
The R and C models are exactly the same versions of the model.
calculate(Cmodel)
## [1] -8089.009
You can simulate from prior.
Cmodel$simulate('lam0')
calculate(Cmodel)
## [1] -7669.788
Specify MCMC.
conf <- configureMCMC(Rmodel,
monitors = c("N", "lam0", "psi", "sigma"))
## ===== Monitors =====
## thin = 1: N, lam0, psi, sigma
## ===== Samplers =====
## RW sampler (402)
## - sigma
## - lam0
## - s[] (400 elements)
## conjugate sampler (1)
## - psi
## binary sampler (200)
## - z[] (200 elements)
Build an executable MCMC.
Rmcmc <- buildMCMC(conf)
Compile in C++.
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)
Run compiled model (do not run th uncompiled model with runMCMC(Rmcmc,100)).
samples <- runMCMC(Cmcmc,100)
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
Explore.
dim(samples)
## [1] 100 4
colnames(samples)
## [1] "N" "lam0" "psi" "sigma"
samplesSummary(samples)
## Mean Median St.Dev. 95%CI_low 95%CI_upp
## N 126.4300000 123.5000000 19.7638866 95.9500000 185.3500000
## lam0 1.0657507 1.0453553 0.1245833 0.7794961 1.2568560
## psi 0.6426589 0.6284448 0.1080914 0.4638732 0.9672855
## sigma 2.2640864 2.2800982 0.1605662 1.7889305 2.4354455
Run compiled model.
samplesList <- runMCMC(Cmcmc,
niter = 10000,
nburnin = 5000,
nchains = 2)
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
samples <- rbind(samplesList[[1]],
samplesList[[2]])
str(samples)
## num [1:10000, 1:4] 31 29 27 23 22 22 19 24 28 29 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:4] "N" "lam0" "psi" "sigma"
Calculate ESS effective sample size (should be at least a thousand). Efficiency measures the tradeoff between ESS and computation time (the ratio). For example, an efficiency of 2 means that for every second you run the algorithm, it produces 2 independent samples.
library(coda)
apply(samples, 2, effectiveSize)
## N lam0 psi sigma
## 61.64045 143.41820 58.91601 44.60710
Produce trace an density plots.
library(basicMCMCplots)
chainsPlot(samplesList,
var = c("N", "sigma", "lam0"))
Display summary stats. Compare to the values used to simulate data, in particuler \(N = 50\), \(\sigma = 5\) and \(\lambda_0 = 0.4\).
summary(samples)
## N lam0 psi sigma
## Min. : 9.00 Min. :0.2070 Min. :0.02896 Min. :2.773
## 1st Qu.: 38.00 1st Qu.:0.4230 1st Qu.:0.19116 1st Qu.:3.965
## Median : 55.00 Median :0.5031 Median :0.27948 Median :4.504
## Mean : 61.08 Mean :0.5151 Mean :0.30760 Mean :4.619
## 3rd Qu.: 80.00 3rd Qu.:0.5903 3rd Qu.:0.40230 3rd Qu.:5.076
## Max. :174.00 Max. :1.3284 Max. :0.87699 Max. :9.349